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We present a method for extracting gravitational waves from numerical spacetimes which gen- 
eralizes and refines one of the standard methods based on the Regge-Wheeler-Zerilli perturbation 
formalism. At the analytical level, this generalization allows a much more general class of slicing 
conditions for the background geometry, and is thus not restricted to Schwarzschild-like coordinates. 
At the numerical level, our approach uses high order multi-block methods, which improve both the 
accuracy of our simulations and of our extraction procedure. In particular, the latter is simplified 
since there is no need for interpolation, and we can afford to extract accurate waves at large radii 
with only little additional computational effort. We then present fully nonlinear three-dimensional 
numerical evolutions of a distorted Schwarzschild black hole in Kerr-Schild coordinates with an odd 
parity perturbation and analyze the improvement we gain from our generalized wave extraction, 
comparing our new method to the standard one. We do so by comparing the extracted waves with 
one-dimensional high resolution solutions of the corresponding generalized Regge- Wheeler equation. 

We find that, even with observers as far out as R = 80 M — which is larger than what is commonly 
used in state-of-the-art simulations — the assumption in the standard method that the background 
is close to having Schwarzschild-like coordinates increases the error in the extracted waveforms 
considerably. Even for our coarsest resolutions, our new method decreases the error by between 
one and two orders of magnitudes. Furthermore, we explicitly see that the errors in the extracted 
waveforms obtained by the standard method do not converge to zero with increasing resolution. 
That is, these errors are dominated by the extraction method itself and not by the accuracy of our 
simulations. We analyze in detail the quasinormal frequencies of the extracted waves, using both 
methods. 

In a general scenario, for example a collision of compact objects, there is no precise definition 
of gravitational radiation at a finite distance, and gravitational wave extraction methods at such 
distances are thus inherently approximate. The results of this paper bring up the possibility that 
different choices in the wave extraction procedure at a fixed and finite distance may result in relative 
differences in the waveforms which are actually larger than the numerical errors in the solution. 



PACS numbers: 04.25.Dm, 04.25. Nx, 04.70.Bw 
I. INTRODUCTION 



One of the goals of numerical solutions of Einstein's 
equations is usually the prediction and analysis of the 
gravitational radiation emitted in some physical process. 
There are many methods for computing, or extracting, 
gravitational waves from a numerical spacetime. They 
can be broadly divided into two groups, depending on 
whether the solution includes null infinity (or a portion 
of it), or whether the computational domain is truncated 
at a hopefully large but finite distance from the source. 
In the first case, gravitational radiation can be defined 
and extracted in an unambiguous, rigorous way (see e.g. 
[l[ and references therein, and Q)- In the second case, 
some approximation has to be made; not only at the 
level of the observer being in the radiation zone, but also 



in the way the "gravitational radiation" is computed in 
terms of the spacetime metric. Due to the additional 
complexity of evolving Einstein's equations all the way 
up to null infinity, currently most simulations actually 
truncate the computational domain by placing an arti- 
ficial outer boundary at a finite distance. This paper 
deals with one particular approach to gravitational wave 
extraction from spacetimes within this second group. 

In general, one expects the differences between the ex- 
act waveforms and those extracted at a finite distance 
to decay as the extraction radius increases. One natural 
question that arises is: for a given extraction method, 
how far away is far enough, so that the errors in the ex- 
tracted waves are dominated by the accuracy of the sim- 
ulations used to obtain the numerical spacetime, and not 
by the extraction mechanism itself? In this paper we ad- 
dress this question in detail in a very particular scenario, 
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but which might shed some light on the general case. 

The main idea of extracting waves at a finite dis- 
tance is to exploit the structure of an asymptotically 
fiat spacetime. One reads off the quantities which are 
needed to compute the gravitational radiation from the 
numerically generated solution. The method which we 
consider in this paper is based on the well-known per- 
turbations of the Schwarzschild spacetime. See e.g. 
(E H S, S, 0, B ©, M, d E El Q for other approaches 
based on the Weyl scalar <3>4. 

One possible approach is to assume that the full metric 
in the region of extraction can be considered as a pertur- 
bation of a flat spacetime, and to read off such pertur- 
bations from the numerical solution. This approach is 
justified by the fact that the leading order of the metric 
at large distances (in an expansion in powers of 1/r) is 
flat. If the waves are extracted at a large but finite dis- 
tance from the source, it makes sense to try to decrease 
the errors of the approximation by further considering 
the next order in the expansion of the metric, which is 
described by the Schwarzschild solution. In doing so, the 
numerical metric is not considered anymore a perturba- 
tion of flat spacetime, but instead of the Schwarzschild 
geometry. One can consider even higher orders in this 
background identification, such as the spin contribution. 
However, an important fact to keep in mind is that all 
these methods should in principle give the same gravita- 
tional radiation as the radius of extraction increases. In 
other words, one should be able to compute the gravita- 
tional radiation through, for example, perturbations of a 
flat spacetime or the Schwarzschild metric, and the ra- 
diation should contain the information about the space- 
time's non-zero mass and — if present — angular momen- 
tum when the observer is at large enough distances. 

If only the first or the first two orders in the asymp- 
totic expansion of the metric are kept when identifying 
this distant "background" geometry, then the framework 
for extracting gravitational radiation is that of perturba- 
tions of flat spacetime or of the Schwarzschild geometry, 
respectively. One can view the former as a subcase of 
the latter, so that from hereon we will just consider per- 
turbations of the Schwarzschild spacetime. In this case, 
perturbations decouple into two separate sectors, which 
differ in the parity of the perturbations (odd or even). 
These two parity sectors are directly related to the real 
and imaginary parts of the Weyl scalar ^4 (see, for ex- 
ample, ref. [l5|). Gauge invariant formalisms for such 



perturbations were developed by Regge and Wheeler 16 1 
in the fifties for the odd-parity sector and by Zcrilli [17 1 
in the seventies for the even-parity sector. 

The idea of using Regge Wheeler Zerilli perturbation 
theory to extract gravitational waves from numerical 
spacetimes is definitely not new. It goe s back to pioneer- 
ing work by Abrahams and Evans |18l . Il9l . |20| (see also 
j2l| ) and it has been used extensively since the birth of 
numerical relativity (see [22j for a review). For example, 
the accuracy of simulations of distorted black holes was 
tested by comparing extracted waveforms against per- 



turbative calculations [23, [2J, [25|, [26J, [27| , and often, also 
technical improvements (such as excision) were tested by 
studying their effects on waveforms [lH, . Recently, 
[30j | reported Zerilli waveforms from unequal mass bi- 
nary black hole inspirals. In hydrodynamical simula- 
tions, gravitational waves are often determined via the 
quadrupolc formula, which usually gives more accurate 
information in these particular situations (unless a black 
hole is present), since the wave amplitude is typically 
very small and thus difficult to detect from the space- 
time metric [HI, [H, [H| . 

In this paper we present a generalization of this ap- 
proach to gravitational wave extraction with two salient 
features. The first is at the level of the perturbation for- 
malism itself: we use a generalization of the standard 
Regge-Wheeler-Zerilli (RWZ) formalism, which is not 
only gauge invariant, but also covariant [HI. l34l. HH, |3(| . 
in the sense that it is independent of the background co- 
ordinates. The standard RWZ formalism is gauge invari- 
ant only in the sense that the background metric is fixed 
to the Schwarzschild geometry in Schwarzschild coordi- 
nates, and the formalism is invariant with respect to in- 
finitesimal, first order changes of coordinates, which keep 
the background coordinates fixed. However, in numerical 
simulations of Einstein's equations, the numerical space- 
time might be close to the Schwarzschild geometry in 
certain situations (say, at large distances), but the met- 
ric does not need to be close to the Schwarzschild metric 
in Schwarzschild coordinates. In fact, when dealing with 
the black hole singularity through black hole excision, 
one uses coordinates that are well defined in a neighbor- 
hood of the horizon, and which are therefore clearly not 
of Schwarzschild type. 

This first salient improvement (the use of a generalized 
formalism) is independent of the details of the numeri- 
cal implementation. The second improvement is tied to 
our particular numerical approach, which uses high or- 
der methods (typically higher than four) for high accu- 
racy, and uses multiple blocks with adapted grids, non- 
trivial topologies, and smooth boundaries. The use of 
high order methods for both the evolution of Einstein's 
equations and for the wave extraction procedure itself, 
combined with the use of shells of "spherical" patches 
or blocks, allows us to extract gravitational waves in a 
simple, fast, and accurate way. In particular, we can 
keep both the angular and radial resolutions fixed and 
place the outer boundaries at large distances, using con- 
siderably less computational resources than what would 
be needed with Cartesian grids, even when using mesh 
refinement. In addition, no interpolation to spheres is 
needed to extract waves on spherical shells. 

For weak perturbations of a Schwarzschild black hole 
we can actually obtain the exact solution by evolving 
the generalized Regge-Wheeler equation. Since this is 
a wave equation in 1 + 1 dimensions, we can solve it 
with almost arbitrarily high accuracy. For all practical 
purposes, we consider it to be an exact solution, against 
which we can compare the extracted waveforms from our 
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three-dimensional evolutions. We evolve weak perturba- 
tions of a Schwarzschild black hole in Kerr-Schild coor- 
dinates, using the fully nonlinear Einstein equations. We 
find that the assumption in the standard method that 
the background is in Schwarzschild-like coordinates in- 
creases the error in the extracted waves (as compared to 
extracting with the correct background) by between one 
and two orders of magnitudes. This is true even with 
observers as far away as 80 M, and even for the coarsest 
resolutions that we use. Furthermore, we explicitly see 
that the errors in the standard method do not converge 
to zero with increasing resolution at any fixed extraction 
radius, while they do with the generalized method. The 
errors only decrease (as 1/r, as we discuss in sect. IIV[) as 
the observer radius is increased. That is, if one does not 
use the correct background coordinates, these errors are 
dominated by the extraction procedure and not by the 
accuracy of the simulations. We compare the quasinor- 
mal frequencies of the waves extracted the above methods 
against the results predicted by perturbation theory. 

The organization of this paper is as follows. In sect. [TT] 
we describe in a self-contained way the generalized per- 
turbation formalism, restricted to the odd-parity sector 
(we will present a similar treatment for the even parity 
sector elsewhere), and our construction of the Regge- 
Whcclcr function from a numerical spacetimc. We also 
use the inverse problem (that is, the generation of a per- 
turbed metric from any given Regge- Wheeler function) 
to construct initial data that automatically satisfies the 
Einstein constraints when linearized around the Schwarz- 
schild spacetime, which does not necessarily need to be 
given in Schwarzschild coordinates. This is the data that 
we later evolve and use in our numerical tests. 

In sect. IIIII we briefly describe our numerical tech- 
niques, out formulation of Einstein's equations, and our 
outer boundary conditions. Finally, we present our nu- 
merical results in sect. IIVI We first show that our ex- 
tracted covariant and gauge invariant Regge- Wheeler 
function coincides very well with the expected one from 
perturbation theory (which we obtain by solving the 1 + 1 
generalized Regge- Wheeler equation) when we use the 
generalized formalism to identify the background cor- 
rectly. After that, we compare our covariant and gauge 
invariant extracted waveforms with those obtained by 
the traditional approach, which assumes that the back- 
ground is either the Minkowski spacetime in Minkowski 
coordinates, or the Schwarzschild spacetime in Schwarz- 
schild coordinates. In sect. |V] we discuss these results in 
the broader context of gravitational wave extraction for 
generic spacetimes. 

For completeness, in Appendix [Al we describe in detail 
our conventions for tensor spherical harmonics decompo- 
sitions. 



II. ODD-PARITY PERTURBATIONS OF 
SCHWARZSCHILD AND WAVE EXTRACTION 

This section summarizes the results of the generalized 
formalism relevant for this paper. We closely follow the 
notation and presentation of ref. [l5j . 

A. The background metric and tensor spherical 
decomposition of the perturbations 

The generalized formalism assumes that the total met- 
ric can be written as 

9™ = 9i«> + 89iw (!) 

where g^ v describes the Schwarzschild geometry and 5g^ v 
is, in some sense, a "small" correction. Further, it is as- 
sumed that the four-dimensional manifold can be decom- 
posed as the product of a two-dimensional manifold M. 
parametrized with coordinates x a (a = 0, 1) and a unit 
2-sphere S 2 with coordinates x A (A = 2,3), such that 
the background Schwarzschild metric takes the form 

ds 2 = g ab {t,r)dx a dx b + f 2 (t,r)g AB dx A dx B . (2) 

Capital Latin indices refer to angular coordinates (9, cj>) 
on S 2 , while lower-case ones refer to the (t,r) coordi- 
nates. Here §ab is the standard metric on the unit 
sphere, g a b denotes the metric tensor on the manifold 
M, and f 2 is a positive function. If one uses an areal 
radius coordinate, then / = r, but we do not make 
such an assumption. Actually, as we discuss below, the 
fact that our formalism is general enough to allow for 
/ = f(t,r) has practical advantages in the wave extrac- 
tion procedure. For simplicity, the metric on the unit 
2-sphere S 2 is assumed to be in standard coordinates: 
fJAB = diag(l,sin 9). Summarizing, we are assuming 
that the background Schwarzschild metric is given in a 
coordinate system in which there is no angular shift, but 
there can be a radial shift. Note that there is no assump- 
tion about the shift in the perturbation. 

From a numerical relativity point of view, it is usually 
convenient to deal with the variables that appear in the 
3+1 split of sp acetime. To this end, we follow the nota- 
tion of ref. [la ] and explicitly expand the components of 
the background Schwarzschild metric as 

ds 2 = (-a 2 +j 2 /3 2 )dt 2 + 2j 2 (3dtdr + j 2 dr 2 (3) 
+ f 2 (d9 2 + sin 2 9 d(f> 2 ) 

where a and (3 = (3 r are the background lapse and ra- 
dial shift vector, respectively, and j 2 = g rr . Since the 
background is spherically symmetric, it is convenient to 
expand the perturbations in spherical harmonics, 

oo t 

§ 9^ = E E ■ ( 4 ) 

£=1 m=-i 
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In the odd-parity sector there is no perturbation for 
I = 0. The dipole term, £ = 1, corresponds to the lin- 
earization of the Kerr metric using the angular momen- 
tum of the spacetime as a parameter. Thus, for gravita- 
tional wave extraction we only need to consider pertur- 
bations with I > 2. These quantities can be parametrizes 
according to 



(5) 
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Using the covariant derivative V a compatible with the 
metric c/ab on the unit sphere S 2 and its associated Lcvi- 
Civita tensor £ab (with non- vanishing components eg^ = 
sin 6 = — £0e ) i wc define Sa = f^VflF (the first index 
in e raised with the inverse of g) and Sab = V(aSb)- 
Here, Y = Y^ m ) are the standard spherical harmonics. 
The quantities Sa and Sab form a basis on S 2 for odd- 
parity vector and symmetric tensor fields, respectively. 
For completeness, we give a detailed and self-consistent 
description of how to use these to decompose vectors and 
tensors into spherical harmonics in Appendix [SJ 

From now on, wc suppress the supcrindices (£, m) and 
the sum over them, since modes belonging to different 
pairs of (£, m) decouple from each other in the perturba- 
tion formalism. 



B. Extraction of the Regge— Wheeler function from 
a given geometry 

To define the background metric, we extract the I = 
component (that is, the spherically symmetric part) of 
the numerical solution 5*°*. This is done by decomposing 
the metric g a b of the two-dimensional manifold M. into 
spherical harmonics. These metric components behave 
like scalars under a rotation of coordinates. Thus, the 
background metric is computed as 



~ 9ab 47 



tot dn , 



9ab 



(G) 



where dfl is the standard area element on S 2 . The func- 



tion / can be computed through / = yA/47r, with 



A 



' g dd d<f> , 



(7) 



where the integration is performed over the extraction 
2-sphcre, and g is the determinant of cjab ■ 

Similarly, we compute the perturbed quantities by ex- 
tracting the I > 2 components of the numerical metric 
5^°', in the way explained in Appendix [XJ 



Once we have obtained the multipolcs &, hi, /12, 7Ti, 112 
defined above in eq. ([S]) and the background quantities 
/, a, 7, (3 defined in cq. J3|), we can find the generalized 
gauge- invariant Regge- Wheeler (RW) function §rw- It 
is given by [l5[ 
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(8) 



where do = dt — f3d r and A = (£— 1)(£ + 2). Notice from 
eq. ([5]) that the only multipole components appearing in 
the RW function Qrw are hi and 7Ti, so that there is no 
need to compute the others. 

Previous approaches to compute waveforms with the 
standard RWZ formalism have typically been consider- 
ably more involved than what we have just described. 
We briefly sketch the standard approach here. Einstein's 
equations are usually solved using Cartesian coordinates 
on a Cartesian grid. The numerically obtained metric 
is first transformed to polar-spherical coordinates. Per- 
forming the multipole decomposition on a given coordi- 
nate sphere requires a numerical integration over that 
sphere, which in turn requires interpolating the metric 
to the spherical surface, which does not coincide with 
the grid points of the Cartesian grid. Integrating over 
the sphere also allows computing the areal radius and 
its radial derivatives. These quantities are then used to 
transform the metric in a second step to its final form 
in "Schwarzchild-like" coordinates. This is done by first 
changing from the coordinate radius to an areal radius 
(which requires the numerically calculated radial deriva- 
tives), and then identifying the (t,r) components of the 
metric in this new coordinate system, which is assumed to 
be a perturbation of the Schwarzschild metric in Schwarz- 
schild coordinates. With all this in place, the waveforms 
are then computed using standard RWZ formulae. 

In our case, the multi-block grid structure naturally 
allows for spherical surfaces. Hence, no interpolation is 
required. The generalized perturbation formalism allows 
us to compute the RW function <&rw without transform- 
ing the metric to Schwarzchild coordinates. In particular, 
the transformation to an areal radial coordinate is not 
required at all. Thus, our extraction procedure amounts 
simply to numerical integrations at a given value of the 
radial coordinate to compute the multipolcs, and then us- 
ing eq. ([5]) to compute the RW function. An additional 
improvement is that our high order accurate derivative 
operators are naturally associated with a high order ac- 
curate discrete norm, leading to an integration procedure 
which has the same accuracy as our derivative operators. 



(Re)construction of the metric from the 
Regge— Wheeler function 



It can be seen (see, for example, ref. 15( for more de- 
tails of what follows) that for any slicing of Schwarzschild 
of the type given in eqs. or ([3]), that we can construct 
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a perturbed four-metric from the RW potential. The per- in eq. ((T|) . become 
turbation coefficients of the linearized metric, as denned 



5g r e 



Jnb 



till 



1 (-f*RW + m'sw + *Rw(Pf - /)) + f — f 2kf ] ^ 
a, \ / j sin a 

1 (-f<s>Rw + m' RW + ®rwW - ft) + fk '~ 2kf 

a V ) I 



sin 6% 



2k 



sin 
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sin 9 
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(10) 

(11) 
(12) 

(13) 

(14) 



sm9Y e (15) 



r 



Here dots and primes denote derivatives with respect to 
time and radius, respectively. It is Y$ — d<pY, Yg = dgY, 
and as before we are skipping the (i,m) superindices. 
7, a, (3, and / are defined in cq. ([3]). It can be seen 
that the function A: is a pure gauge term and completely 



arbitrary; in particular, we can make it vanish (resulting 
in the so called Regge-Wheeler gauge) through a first 
order coordinate transformation. 

The generalized RW equation is 



*flff = ci&rw + c 2 ® RW + cz^rw + C4&RW - ct 2 V<& RW 
with the coefficients Cj and the potential V given by 

ci = 2/3 

(« 2 -7 2 /? 2 ) 



C-2 



C3 = 



(7a — 7/3a' + afij' — cry + ^af}') 



c 4 = — — -Y(ia - cry + Y(3 2 a' - 2^aPf3' + 7 d cr3 + T^Pi + 7crV - p a Pl' 
7 J a 



V = —r 
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P 



1) - 



6M 



(16) 

(17) 
(18) 

(19) 

(20) 

(21) 



When the background metric is Schwarzschild in 
Schwarzschild coordinates, this generalized RW equation 
coincides of course with the standard equation. Below, 
in sect. II VI we use high-resolution solutions of this gener- 
alized 1 + 1 equation as "exact" solutions , against which 
we compare the extracted RW function from our three- 
dimensional distorted black hole simulations. 



III. FORMULATION OF THE EQUATIONS, 
BOUNDARY CONDITIONS, INITIAL DATA, 
AND NUMERICAL METHODS 

A. Evolution equations 

The numerical simulations shown below were per- 
formed by evolving a first order symmetric hyperbolic 
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reduction of the Generalized Harmonic formalism, as con- 
structed in ref. [3^ |, In this formulation, the coordinates 
x M are chosen to satisfy the (generalized) harmonic con- 
dition 1 m 

V ff V n s W = H^\t,x l ), (22) 

where the gauge source functions H^{t,x l ) are freely 
specifiable functions of the spacetime coordinates, and 
V M is the covariant derivative associated with g^. Here 
we omit the label "tot" from the metric ([I]) for the sake 
of simplicity. The reduction from second to first order is 
achieved by introducing the first derivatives of the metric 
g^ v as independent quantities. Following ref. [37j | . we 
introduce the quantities 

Qp, v = -n a d a g^ (23) 
D HIV = dig^, (24) 

where = —aV^t is the (future directed) timelike unit 
normal vector to the hypersurfacc t — x° = const. Thus, 
the evolution equations for Q^u are given by the General- 
ized Harmonic formalism, while the evolution equations 
for are obtained by applying a time derivative to 

their definition (|24| and commuting the temporal and 
spatial derivatives. Finally, the metric g^ is evolved us- 
ing the definition of Q^ v , eq. (|23[) . In addition, in the 
spirit of refs. [39|, |4(J, the constraints of this system are 
added to the evolution equations in such a way that the 
physical solutions (i.e., those satisfying the constraints) 
are an attractor in certain spacetimes. In those situa- 
tions, small constraint violations will be damped during 
the evolution. The whole construction of this formulation 
of the equations is described in detail in [33| ■ 

The standard 3 + 1 components of the metric (i.e., the 
lapse function a, the shift j3 z , and the intrinsic metric 



7y) can be obtained via the relations 

a 2 = -l/g u (25) 
ft = 9ti (26) 
lij = 9ij ■ (27) 

The extrinsic curvature is defined in terms of the intrinsic 
metric as 

Kii = ^{d t ~Cp) ll3 . (28) 

It can be recovered from the fields Q^, via 

-2aK ij = a(),, ■ I)„, ■ I),,, (29) 
-1 km f3 m (D ijk +D jik ). 

1 In this subsection we use the Latin indices i,j,k,. . . to denote 



three-dimensional spatial quantities, while Greek indices con- 
tinue to represent the four-dimensional ones. 



In our simulations below we also monitor the Hamilto- 
nian and momentum ADM constraints, namely, 

H = l^R-K^ + itrKf) (30) 
M t = V fe (If* - SftrK) , (31) 

where ( 3 'R is the Ricci scalar associated to the three- 
dimensional space-like metric 7^ . 

We impose maximally dissipative boundary conditions 
at the outer boundary. While these conditions guaran- 
tee well-poscdness of the associated initial value prob- 
lem, and thus numerical stability with our particular dis- 
cretization, they are physically incorrect in the sense that 
they do not include back-scattered radiation from outside 
the simulation domain. For that reason, in the simula- 
tions shown below we place the outer boundary at large 
enough distances so that our extracted waves are causally 
disconnected from boundary effects. 



B. Multi-block approach 

We use multi-block (also called multi-patch or multi- 
domain) methods for our numerical calculations. These 
have several advantages over single-domain methods: 

Smooth boundaries. They provide smooth outer and 
inner boundaries, which is in general required [4~0 | for a 
well-posed initial boundary value problem. 

Constant resolution. They allow us to use constant 
radial and angular resolutions. This is not possible with 
mesh refinement methods. The way in which mesh re- 
finement is typically used leads to a decreasing radial 
resolution, which makes it difficult to extract accurate 
gravitational wave information in the wave zone of a bi- 
nary black hole system. (See sect.lVl where we list typical 
wave extraction radii and resolutions.) 

No time-stepping restrictions. They do not lead to a 
deterioration of the CFL factor for co-rotating coordi- 
nates. (For example, [1^ reports that the CFL factor 
had to be reduced on the outermost refinement levels. 
The same was done later in [H|].) 

Adapted to symmetries. They can be adapted to the 
symmetries of a system. Obviously, adapted coordinates 
can reduce the discretization error significantly. In our 
case, we use on each block one radial and two angular 
coordinates to model the geometry of a single black hole. 
For binary black hole systems, one can use blocks that are 
roughly spherical near the individual holes and far away 
in the wave zone, with a transition region in between. 
Fig. 5 in [3| shows a possible multi-block system for 
this. 

No coordinate singularities. They have no coordinate 
singularities. Spherical or cylindrical coordinates have 
singularities on the z axis which may cause problems. 
An alternative approach which avoids these singularities 
would be to use a pseudo-spectral decomposition into 
spherical harmonics. This was used in [45j to evolve 
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which may be different. These penalty terms ensure con- 
tinuity between the two blocks in the continuum limit and 
numerical stability in the semidiscrete case if the relevant 
parameters are appropriately chosen. The quantities a 1 
and cr r depend on the coefficients of the differencing op- 
erators that are used on the two blocks. 2 The parameters 
S l and S r determine how much (if any) dissipation is in- 
troduced across the block boundary. To ensure stability, 
they must be chosen in a very specific way depending on 
the characteristic speeds of the evolution system. 

We have implemented this in the Cactus framework 
HI, HH , using the Carpet driver [64], [fill and the Cactus- 
Einstein toolkit 1661] . 



Figure 1: The equatorial plane of an example six-block geom- 
etry, cutting through four blocks. Note that the blocks do not 
overlap. All six blocks are made up identically. The outer and 
inner boundaries are smooth spheres. The outer boundary in 
our typical simulations is actually located much further out 
than shown here. 



scalar fields on a Kerr background, in 46j, [47J to evolve 
Einstein's, and in [H, [49|, Ho, [EH, [E2, [53 1 to set up initial 
data for various black hole and neutron star configura- 
tions. 

Of course, using multiple blocks adds to the complex- 
ity of an implementation. However, the properties of 
multi-block systems for hyperbolic equations are by now 
well understood, and we describe our particular approach 
in [54| and [55|, and in some detail in [44[. In this 
particular paper we use a six-block system to discretize 
the geometry of a single black hole, which is depicted 
in fig. [TJ We use the same tensor basis on each block, 
namely a global three-dimensional Cartesian coordinate 
system. We found that this greatly simplifies the inter- 
block boundary conditions, since all components of tenso- 
rial quantities are then scalars with respect to the block- 
local coordinate systems. 

We use the penalty method to enforce the inter-block 
boundary conditions. The penalty method for finite dif- 
ferences is described in [EH, |57l IBH . and we describe our 
approach and notation in [44l. |54|. [55| . In short, the 
penalty method works as follows. The individual blocks 
do not overlap, but they have their boundary points in 
common. The evolution equations are first discretized on 
each block independently using one-sided derivatives near 
the block boundaries. Then a correction term is added 
to the right hand side of the time derivative of each char- 
acteristic variable at the boundary points, penalising the 
difference between the left and right eigenmode values u l 
and u r on the boundary points: 

<?' 

cV -► cV + -r-r (u r - u l ) (32) 
h l a L 

d t u r -► d t u r + {u l - u r ) . (33) 
h r a r K ' x ' 

Here h l and h r are the grid spacings on the two blocks, 



C. Initial data 

If the RW function $rw satisfies the RW equation 
(fl6|) . then the perturbed metric constructed in sect. Ill Cl 
satisfies the linearized Einstein equations. Furthermore, 
it can be explicitly shown that this metric initially sat- 
isfies the linearized constraints around the Schwarzschild 
geometry for any initial values <&Rw(t = 0, r) and 
§Rw(t = 0, r). 3 We take advantage of this property and 
construct initial data in a simple way as a test our new 
wave extraction method. For our simulations below, we 
use Kerr-Schild coordinates for the Schwarzschild back- 
ground, and for the distortion we set £ = 2, m = 0, and 
choose 



$Rw(t = 0,r) = 0, 

$ RW (t = 0,r) = Ae {r - ro)2/a2 



(34) 



with parameters tq and a. This corresponds to a Gaus- 
sian pulse of width a centered at r = tq. 

If we assume that we can Taylor-expand (a suitable 
norm of) the discrete non-linear constraints in terms of 
the perturbation amplitude A for any fixed gridspacing 
h, we have 



C(A, h) 



C(A, h)\ A=Q 
dC{A, h) 



(35) 



.4 
2 



OA 



d 2 C(A,h) 



A=0 



OA 2 



+ 0(A 3 ). 



A=0 



Since in the continuum the linearized constraints are sat- 
isfied, the first two terms in the above expansion van- 
ish for h — > 0, but otherwise are of the order of the 



2 To be exact, a 1 and a r depend on the coefficients of the discrete 
norms that are used in the blocks, but the differencing operators 
and the norms are usually chosen together to satisfy summation 
by parts. See [59I . kill |61| , and especially |55l . 

3 When constructing initial data for the 3 + 1 quantities, one also 
needs to take time derivatives of the four-metric; where second 
time derivatives of Qrw appear, we use the RW equation to 
trade these for space derivatives. 
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truncation error. For small enough A the first term 
(that is, the background contribution) dominates, and 
the term C(A, h) appears to be independent of A. For 
large enough A, on the other hand, the quadratic term 
in the expansion given by eq. (|35[) will dominate. 

Fig. [5] presents numerical evidence that this expected 
behavior is indeed the case. We set up numerical data 
according to eq. (|34| . with perturbation amplitudes A 
between 10 -6 and 10 _1 . The radial domain extent is 
1.8 < r < 7.8, the perturbation is centered around 
r = 4.8 M and has a width of a = 1.0 M. We then 
compute the discrete Hamiltonian and momentum con- 
straints H and M l for these initial data sets, using the 
same (high) resolution, namely 109 x 109 grid points on 
each block in the angular direction and 406 points in the 
radial direction, corresponding to Ar ps 0.0148 M. Due 
to the symmetry of our six-block structure and the ax- 
isymmetry of the initial data, two components of the dis- 
crete momentum constraints coincide, A4 X = A4 y , and 
we therefore do not show the latter. The behavior of 
the constraints in the L2 and the L^, (not shown in the 
figure) norms agrees with eq. (|35[) : for small amplitudes 
A, the discrete constraints at a fixed resolution appear 
to be independent of A, while for large amplitudes they 
show the expected quadratic dependence on A. We also 
show that the discrete constraint violations of our ini- 
tial data sets have the expected dependence on resolu- 
tion. For small amplitudes and coarse resolutions, the 
contribution of the quadratic term in eq. (]35[) is suffi- 
ciently small, so that the constraints seem to converge 
towards zero. However, for any given amplitude A a fine 
enough resolution h reveals that the convergence is actu- 
ally towards a small but non-zero value, determined by 
the quadratic term in the expansion eq. (|35p . This be- 
havior is shown in fig. [3] As an illustration we show there 
a convergence test for H by comparing initial data for dif- 
ferent resolutions. The highest resolutions are identical 
to those used in fig. [5] The other four resolutions shown 
are 73x73x271, 49x49x181, 25x25x91, and 17x17x61 
grid points per block, corresponding to Ar « 0.0222 M, 
0.0333 M, 0.0667 M, and 0.1 M, respectively. 



IV. NUMERICAL STUDIES 
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Figure 2: Discrete constraint violations for various perturba- 
tion amplitudes A at a fixed (high) resolution. We show the 
L2 norm for the Hamiltonian constraint and for two compo- 
nents (a; and z) of the momentum constraint (which turn out 
to be very close to each other, as the plot shows). The nu- 
merical resolution is 109 x 109 grid points per block in the 
angular directions and Ar ~ 0.0148 in the radial direction. 
The behavior is as expected and as described in the body of 
the paper: for sufficiently small amplitudes, the background 
contribution dominates the discretization error in the con- 
straints, which then appear to be independent of A. For large 
enough amplitudes, the constraint violation has a quadratic 
dependence on A (with an exponent of 2.01T0.01 for the reso- 
lution shown in this figure), since for our initial data only the 
linearized constraints (around Schwarzschild) are satisfied. 
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A. Description of the simulations 

We use the Z?8-4 operator constructed in [5{|, a sum- 
mation by parts operator [oTI . [q8| which is eighth order 
accurate in the interior and fourth order accurate at the 
boundaries, optimized to minimize its spectral radius and 
boundary trunctation errors. Fifth order global conver- 
gence is expected [H, IzO| - 

We integrate in time with a 
fourth order Runge-Kutta integrator with adaptive time 
stepping as described in [71 1. 

In order to test both the long term stability and the 
convergence of our code, we first evolve a Kerr black hole 
in Kerr-Schild coordinates with spin j = 0.5. Fig. Q] 



Figure 3: L2 norm of the Hamiltonian constraint for different 
amplitudes A of the perturbation and for different resolutions 
h. The coarsest resolution uses 17 x 17 points per block in 
the angular directions and Ar = 0.1 M in the radial direction. 
We increase the resolution in all directions, up to 109 x 109 
points in the angular directions and Ar ~ 0.0148 M in the 
radial direction. Since only the linearized constraints are sat- 
isfied, the non-linear constraints do not converge to zero. For 
sufficiently large perturbation amplitudes and for sufficiently 
fine resolutions, the non-linear effects become visible, and the 
constraint violations converge to a constant value which de- 
pends on the amplitude A. As shown in the previous figure, 
this dependence is quadratic, as expected. 
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One of the goals of the analysis that follows is to study 
the effect of the choice of the background metric on the 
accuracy of the waveforms. Since for this scattering prob- 
lem solutions in closed form are not known, we compare 
the waves which we extract from our three-dimensional 
simulations to results obtained with an independent 
fourth order accurate one-dimensional code which solves 
the Reggc- Wheeler equation ([To]) . These ID results were 
obtained with a resolution of Ar = 0.0125 M. The rela- 
tive difference in this Regge- Wheeler function to a result 
from twice this resolution lies roughly between roundoff 
error and 10 -7 , which is far below the numerical errors 
that we expect from our 3D simulations. Therefore, we 
consider these ID results in the following to be exact for 
all practical purposes. 




4.8 1 1 1 1 1 

500 1000 1500 2000 

time / M 

Figure 4: L2 norm (top panel) and convergence factor (bot- 
tom panel) for the Hamiltonian constraint for evolutions of a 
Kerr black hole with spin j — 0.5. The coarse resolution cor- 
responds to 16 x 16 points per block in the angular directions 
and Ar = 0.2 M in the radial direction. The fine resolution 
a factor of 1.5 higher in all directions. We see fifth order 
convergence, as expected for the difference operators used. 



shows the L2 norm of the Hamiltonian constraint vs. time 
for two different resolutions. The radial domain extent is 
1.8 M < r < 11.8 M. The coarse resolution corresponds 
to Ar = 0.2 M and 16x16 points per block in the angular 
directions, and the fine resolution increases the number 
of points in all directions by a factor of 1.5. We see 
approximate fifth order convergence, as expected. 

In the simulations discussed below, we place our in- 
ner boundary at r = 1.8 M and our outer boundary at 
r = 251.8 M. This allows for observer locations up to 
r = 80 M, which are still causally disconnected from the 
outer boundary for times long enough to follow the ring- 
down, namely up to t = 280 M. We set up initial data 
according to eq. (JMJ) with A = 0.01, a = 1.0 M, and 
ro = 20 M, where M is the mass of the black hole when 
the perturbation is switched off. Our coarse resolution 
uses 16 x 16 points per block in the angular directions 
and 1251 points in the radial direction, corresponding to 
Ar = 0.2 M. Our fine resolution uses 1.4 times as many 
grid points in all directions. 



B. The standard and generalized RW approaches: 
numerical comparisons 

We now analyze the results of evolving distorted black 
holes as described above and extracting gravitational 
waves with different methods. 

Fig. [5] shows Regge- Wheeler functions for observers 
at r = 20 M, 40 M, and 80 M, extracted with both our 
generalized approach and the standard one. The data 
have been scaled by a factor of 100 to normalize to an 
initial data amplitude A = 1 in eq. (f3~4")) . Recall that 
we used weak waves of amplitude A = 0.01 for these 
simulations to avoid non- linear effects, and to be able to 
compare with the exact solution, which is only known in 
the linear regime. 

Five waves are shown in fig. [5] for each observer loca- 
tion. Apart from the exact solution, we show two results 
obtained from our generalized approach, which coincide 
with each other in the continuum limit. They differ in 
how the background metric is computed: in one case we 
use the exact expressions for the Kerr-Schild background, 
and in the other case these coefficients were numerically 
calculated by extracting the £ = part of the metric, as 
explained in sect. Ill Bl 

Finally, two waveforms were extracted using the stan- 
dard approach with two different assumptions for the 
background, as found in the literature: a Minkowski 
spacetime in Minkowski coordinates, and a Schwarzschild 
spacetime in Schwarzschild coordinates. We want to 
highlight an interesting feature which can easily be seen 
in cq. j8|). For any observer location, the waves extracted 
with these two background should differ only by a factor 
which is constant in time: 

*w = k*rw, (36) 

where k 2 = gff 1 is radial component of the Schwarzschild 
metric in Schwarzschild coordinates. Such a simple rela- 
tionship is a direct consequence of the vanishing radial 
shift for these backgrounds. We confirmed this expected 
behavior numerically with high accuracy: at all times and 
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Figure 5: Extracted waveforms for observers at 20 M, 40 M, 
and 80 M. Shown is the Regge- Wheeler function obtained 
from the standard RW approach and our generalized one. For 
the former we assumed both a Minkowski background and a 
Schwarzschild background in Schwarzschild coordinates, la- 
beled as RW Mm and RW Sch, respectively. For the gener- 
alized approach we show the results for two cases, in which 
the background metric is dynamically computed from the nu- 
merical solution (Generalized RW I), and where we prescribe 
it analytically (Generalized RW II). Also shown is the exact 
waveform. These simulations were performed with a resolu- 
tion of 16 x 16 grid points in the angular directions on each 
block and Ar = 0.2 M in the radial direction. See the main 
text for more details. 



for all observers we recover this expected ratio between 
the two waves to double precision roundoff error. 

Figure O suggests that, as expected, the differences be- 
tween waves extracted with different methods decrease as 
the extraction radius increases. At r — 80 M, the curves 
show excellent agreement in the L e norm 4 . For a more 
thorough comparison, we look at the differences between 
the extracted waves and the exact solution in fig. [6] For 
consistency with fig. [5l we also scaled the errors relative 
to the initial amplitude of the perturbation. 

Perhaps the most notable feature in fig. [5] is that the 
differences between the waves obtained from generalized 
approach cither with a numerically obtained background 
metric or with the exact (Kerr-Schild) background met- 
ric are smaller than the difference to the exact solution. 
For all practical purposes we can therefore consider them 
identical to each other, and for the rest of the paper we 
leave the latter out of the discussion. 

Fig. [5] also shows that the standard approach — with 
either a Minkowski or Schwarzschild background — leads 
to errors which are considerably larger than the errors 
in our generalized approach, even for an observer at r = 
80 M. For the specific resolution that we used for fig. [BJ 
the errors at r = 20 M with the standard method are 
roughly three orders of magnitude larger than the errors 
with the generalized method. For r = AO M and 80 M, 
the ratio of the errors is of the order 10 3 to 10 1 and 10 2 
to 10°, respectively. 

The previous discussion only analyzes the errors in- 
troduced by the standard method at a fixed resolution. 
Next we discuss the dependence of these results on the 
resolution. It turns out that the difference between the 
different methods is even more striking for higher resolu- 
tions. By construction, the generalized wave extraction 
method should give the exact waveform in the contin- 
uum. At the discrete level, its associated errors should 
converge away with increasing resolution. Fig. [7] shows 
that this is actually the case. On the other hand, the 
errors in the standard approach do not converge to zero, 
as shown in fig. In other words, the accuracy of the 
extracted waves with the standard method is dominated 
by the extraction procedure and not by the numerical 
resolution. 

Fig-El as well as the second panel of fig.[6]show another 
interesting feature. Contrary to expectation, assuming 
Schwarzschild-like coordinates instead of a Minkowski 
background does not necessarily lead to smaller errors in 
the waveforms. For example, for an observer at r = 40M 
and during the time interval of about 25M < t < 50M, 
the errors are actually up to one order of magnitude 
larger for the Schwarzschild-like coordinates. However, 
as can be seen from fig. |6l this feature depends on the 
observer location. We assume that this feature is only a 
coincidence. 



4 Also denoted by -I/ e yeball 
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Figure 7: Shown is a convergence test for the simulations 
presented in the previous two figures. The plots labeled with 
"low res" coincide with the ones shown in the previous figures, 
while the plots labeled with "high res" correspond to 1.4 times 
that resolution. The error in the generalized wave extraction 
method, which by design gives the correct waveform in the 
continuum for these simulations, converges towards zero as 
expected. On the other hand, the errors in the standard wave 
extraction method are almost unaffected by the increased res- 
olution. This indicates that these errors are dominated by the 
extraction method itself, not by the numerical truncation er- 
ror. These results correspond to an observer at 40 M, but they 
look similar for the other extraction radii that we consider in 
this paper. 



Figure 6: Errors for the waveforms shown in fig. [5] 



The plateau in the errors seen in the last 100 M to 
200 M in fig. [6] is due to an offset in the waveform. Wc 
found that, once the wave function decays to a small 
enough amplitude, it no longer oscillates around zero, 
but instead oscillates around a certain offset. This can 
be seen more clearly from the top panel in fig. [5] This 
offset is present for both the standard and the general- 
ized extraction methods; however, there are important 
differences. The first is that the offset for the generalized 
extraction the offset converges to zero with increasing 
resolution, unlike for the standard method. The other is 



that the offset for generalized method is orders of mag- 
nitude smaller than for the standard method. As we 
will discuss in the next subsection, that has direct con- 
sequences when attempting to extract quasinormal fre- 
quencies. This offset is reminiscent of the one that is 
present in RWZ waveforms when there is spin [72l . |73| . 

The oscillatory feature of the wave can be followed for 
a longer time if the offset is subtracted from the waveform 
by hand, that is, if the wave is shifted along the vertical 
axis so that it oscillates around zero at late times. We 
do so by fitting the data to an exponentially decaying 
wave with an offset. (Details about the fit are given in 
the following subsection) The actual values that we de- 
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tcrmined for the offset are given in tabic HI As expected, 
the offset is decreasing with increasing radius for both 
standard RW wave extraction methods. This offset is 
mainly a result of the wrong assumption about the back- 
ground metric, not of numerical error. There is no such 
clear dependence on the radius when using the new gen- 
eralized extraction. Here the offset originates solely from 
truncation error, and converges to zero with increasing 
resolution. This behavior can also be seen in fig. [7] 



Table I: Values of the offset for different wave extraction meth- 
ods and observers at 20M, 40M and 80M. 
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In fig. [5] we show the difference between the waveforms 
shifted by different offset values and the exact solution, 
for the same observers as before. As can be seen from the 
figure, our qualitative statements about the accuracies of 
the different wave extraction methods remain unchanged, 
if you consider the time span during which the amplitude 
of the wave is significant. 5 We conclude that the main 
errors in fig. [6] are not caused by an overall offset in the 
whole waveform. 



C. Quasinormal frequencies 

We now turn our attention to extracting quasinor- 
mal frequencies from the waveforms just discussed. The 
primary goal is to find out whether these frequencies 
are affected by the choice of a specific wave extraction 
method, which may have some presumably small but 
non-vanishing systematic error for any finite extraction 
radius, and if so, by how much. We used data from the 
lower resolution run that we already analyzed in the pre- 
vious section. The accuracy of the frequency does not 
change significantly if we use the higher resolution run 
instead. 

The angular part of the initial data is a pure £ = 2, 
m = mode. Since the background has no angular 



5 Of course, because we subtract the offset by hand to decrease 
the errors at late times, we can naturally follow the oscillatory 
part of the wave for longer times before. 
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Figure 8: Shown are the same quantities as in fig. [6l except 
that an offset is subtracted from each waveform before calcu- 
lating the errors. See the main text for details. 



momentum, there is no mode-mode coupling at the lin- 
ear level, while nonlinear coupling can be neglected for 
the current study, because we only evolve weak pertur- 
bations. Therefore the only dominant multipolc mode 
present in the data at all times should be the one injected 
initially. At the numerical level, I = 4 modes can be gen- 
erated by our six-block grid structure. However, in |74| 
it was found that in the absence of angular momentum, 
these modes not only converge to zero with resolution, 
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but arc also very small for the resolutions considered in 
this paper. In the above reference and in [75j it was also 
shown that overtones are not significantly excited unless 
the black hole is very rapidly rotating. Based on all this, 
we only fit for a single I = 2, m = mode: 

* f *r = Asm(u r t + x) e iWl( *- to) - £ (37) 

where A is the excitation amplitude, uj = ui r + iu>i is the 
complex quasinormal mode frequency, x i s a phase shift, 
£ is the offset and to is the starting time of the quasi- 
normal ringing regime. The latter is not unambigously 
defined (the so called "time-shift problem"), and as a 
consequence neither are the amplitudes of quasinormal 
modes. In ref. [74| it was proposed to minimize the un- 
certainties due to this time-shift problem by looking at 
carefully chosen relative amplitudes (see that reference 
for details). In order to fit numerical data to eq. ([37]) . 
we fix to to an educated guess 6 and then fit for u, A, x, 
and £. Any difference in to is absorbed in A (in which wc 
are not interested at this point) and does not change the 
other extracted parameters. We find the time-window 
of optimal fitting by looking for a local minimum in the 
relative residual between the original waveform and its 
fit. In ref. (74j it was found that such a local minimum 
is usually quite sharp and therefore gives a good criteria 
for choosing the window of time where the quasinormal 
ringing dominates. Similarly, we use the uncertainties in 
this minimum to quantify the errors in the parameters 
obtained in the fit. More details about the fitting proce- 
dure that we use to extract quasinormal parameters are 
given in ref. (74|. 

In the previous subsection wc discussed the presence 
of an offset in the extracted waves with the standard 
method. If such an offset is not taken into account when 
fitting for the quasinormal frequencies (i.e., for a fixed 
£ = 0), eq. ([57)) does not represent the behavior of the 
numerical data well enough, and no reasonable results 
can be obtained from the fit. This is especially the case 
at medium to late time intervals when the amplitude be- 
comes smaller than the offset, so that the wave does not 
cross zero any more. When one tries to fit for these cases, 
the obtained frequency has no relation at all to the cor- 
rect QNM frequency. For example, at r = 20M the offset 
in the waves obtained from the standard RW wave extrac- 
tion is of order 10~ 2 for both a Minkowski background 
and for Schwarzschild-like coordinates. Without taking 
the offset into account, the value of ui r that the fit de- 
termines lies between 10 -14 and 10 -4 , and is of order 
10 -3 to 10 -6 (compare to tablcfll]). In contrast, the offset 
resulting from the generalized RW wave extraction is of 
order 10 -5 for this resolution. This is small enough that 
the problems described above do not play a noticeable 
role. 



For example, taking into account where the initial data and ob- 
server are located, and assuming a propagation speed of one 



Table |TT] shows the complex quasinormal frequencies 
that we obtained from the generalized and from the stan- 
dard RW methods. As mentioned above and discussed 
in detail in ref. [74j|, the error bars are estimated from 
changes in the frequency when changing the time inter- 
val that we use for the fit of the waveform. We assume 
that the predicted frequency from perturbation theory 
for the fundamental £ = 2, m = mode is exact because 
we use a small amplitude for our perturbation. This fre- 
quency is known to be w exact = 0.37367 — 0.08896i (see 
for example [76|). The frequency obtained from the new 
generalized wave extraction is consistent with this exact 
value within the accuracy to which we can obtain these 
numbers from the fit itself. For the standard wave ex- 
traction method, we only find agreement to three signif- 
icant digits in the real part, but better agreement with 
the exact value in the imaginary part of the waveform. 
Note that, since the waveforms only differ by a constant 
factor (see subsec. I IV Bp , the frequencies obtained with 
a Minkowski and a Schwarzschild background agree to 
roundoff error. The reason for the lower accuracy in the 
real part of uj might be due to the fact that the wave- 
forms are slightly distorted due to the wrong assumption 
for the background metric. This causes a larger residual 
between the data and the fit — it is about a factor of two 
larger than with the generalized wave extraction — and 
some degradation in how accurately certain fitting un- 
knowns like ui can be determined. That may also explain 
the larger relative error for the waves extracted with the 
standard RW wave method, which is shown in the right 
column of the same table. There the relative error is 
defined as \(ui — u) cxact )/uj 

exact | ■ 



V. FINAL COMMENTS 

When considering methods for extracting gravitational 
waves from numerical spacctimcs at a finite distance, one 
question of direct interest is: How sensitive is the ac- 
curacy of the extracted waveforms to both the extrac- 
tion method and observer location? In particular, how 
far away is "far enough" when extracting gravitational 
waves? 

It is in general not easy to pose such a question in a 
precise way, since in order to quantify this one needs an 
exact waveform to compare with. This exact waveform 
is in principle only well defined at future null infinity. 
However, there are some particular scenarios of interest 
where the concept of "exact waveforms" at a finite dis- 
tance can be given a well defined and precise sense. That 
is the case, for example, for perturbations of Kerr black 
holes (actually, of Petrov type D spacetimes): the Weyl 
scalar "J 4 is defined everywhere in an essentially gauge 
and tetrad invariant way [73]. Similarly, for perturba- 
tions of Schwarzschild black holes, the Regge- Wheeler 
and Zerilli functions arc defined in a gauge invariant way 
everywhere as well. In fact, there is a one-to-one mapping 
between these functions and ^4; see e.g. (l5| . 
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Table II: Quasinormal frequencies of the I = 2, m = mode as measured by an observer at r — 20 M. Results are given for 
waveforms resulting from the different extraction methods we use. The predicted frequency from perturbation theory, which 
we assume to be exact because our perturbation amplitude is small, is LJ CX act = 0.37367 — 0.08896i [76|]. The uncertainties in 
the extracted frequencies originate from variations in them depending on which interval of the waveform is used for the fit. 
The relative error is defined as | (uj — inexact ) / inexact | . 



Extraction Method 




relative error 


Generalized RW 
RW Min 
RW Sch 


0.3736 - 0.0890i ± (3 + 3i) x 10" 4 
0.3733 - 0.0889i ± (3 + 3i) x 10~ 4 
0.3733 - 0.0889i ± (3 + 3i) x 10" 4 


1.9 x 10" 4 +4.5 x 10~ 4 i 
9.9 x 10" 4 +6.7 x 10~ 4 i 
9.9 x 10" 4 + 6.7 x 10~ 4 i 



Therefore the above question can be posed in a setting 
that might not be the most general one, but it is one in 
which a precise, quantitative answer can be found. The 
concrete setting that we chose for the current study is 
that of weak perturbations of Schwarzschild black holes. 
Furthermore, in this paper we restricted our treatment to 
odd parity perturbations (the even parity sector will be 
presented elsewhere) . One of the standard methods that 
has been widely used for extracting gravitational waves 
from such spacetimes is through the standard Regge- 
Wheeler-Zerilli perturbation formalism. This formalism 
provides a gauge invariant treatment of perturbations 
of a background geometry defined by the Schwarzschild 
spacetime in Schwarzschild coordinates. That is, the 
formalism is invariant with respect to linear coordinate 
transformations which leave the background coordinates 
fixed. If one extracts waves using this formalism from 
a perturbation of Schwarzschild in, say, Kerr-Schild co- 
ordinates, the extracted waves at a finite distance will 
not be correct, even with if extracted with infinite nu- 
merical precision. There is a systematic error in such 
an extraction, due to the incorrect identification of the 
background coordinates. Of course, one expects this er- 
ror to decrease as the extraction radius increases. In 
the spirit of the above discussion, the question that we 
asked ourselves was: how far away must the observer be, 
so that the difference between the exact waveform and 
the extracted one is negligible, if the extraction method 
has a systematic error? For this study we chose a very 
specific interpretation of "negligible systematic errors in 
the waveforms" , namely, that they are smaller than or 
comparable to the errors in these waveforms due to the 
numerical discretization. 

In order to provide a quantitative answer to this ques- 
tion, in this paper we first proceeded to generalize the 
standard Regge- Wheeler extraction approach by using a 
perturbation formalism that allows for quite general slic- 
ing conditions for the Schwarzschild background. With 
this generalization, if one calculated with infinite reso- 
lution, one would extract the exact waveforms for any 
(not necessarily large) finite extraction radius. This holds 
even if the Schwarzschild background is, for example, 
given in a time dependent slicing, or one in which the 
coordinates in a neighborhood of the horizon are well 



defined, as is usually the case in numerical black hole 
evolutions. 

After summarizing the basics of the generalized formal- 
ism, we described our numerical implementation of the 
generalized extraction mechanism and our way of solv- 
ing Einstein's equations. For the latter we used multiple 
blocks and high order methods, both of which present 
several advantages. Of particular interest to this pa- 
per is that, due to the adaptivity that multiple blocks 
provide, the outer boundary can be placed at large dis- 
tances, with much smaller computational costs than with 
Cartesian grids and mesh refinement. We made use of 
this specific advantage and performed three-dimensional 
non-linear simulations of weakly distorted Schwarzschild 
black holes, from which we extracted waves at dis- 
tances larger than most current state of-thc-art three- 
dimensional simulations of Einstein's equations. 

Then we studied the dependence of the extracted wave- 
forms on the extraction method. More precisely, we com- 
pared the standard RW method with our generalized one. 
We found that, even for the coarsest resolutions that 
we used, the errors in the waveforms from the standard 
method were dominated by the extraction procedure and 
not by the numerical accuracy of the spacetime met- 
ric. Furthermore, by increasing the resolution we could 
explicitly demonstrate that the errors in the standard 
method do not approach zero, while they do with the 
generalized one. While this is obviously the expected be- 
havior on analytical grounds, we emphasize that we could 
explicitly sec these differences even with an extraction ra- 
dius which is significantly larger than those typically used 
in current state-of-the-art three-dimensional simulations. 

What is not clear, however, is whether the wave zone 
resolution currently used by mesh refinement codes is suf- 
ficient to sec the differences that we have demonstrated 
in this paper. For example, the spatial resolutions in the 
wave zone of current binary black hole inspiral and co- 
alescence simulations are usually much coarser than the 
resolutions that we used above. Some radial resolutions 
h in the wave zone of binary black hole system simu- 
lations are: [7§| h = 0.5 M . [3f]j h = 0.5 M (but the 
extraction is performed very close in at R = 16M) [7^| 
h = 0.75 M (but h = 1.5 M for calculating the radiated 
angular momentum J), [13] h = 0.85 M, [sij] h = 0.87 M 
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[| h = 0.82 M, [| h = 0.56 M, H ft = 0.56 M. Some 
of these codes are 4th order accurate, but many have at 
least certain components that are only 2nd order accu- 
rate. 7 

One of the interesting features of the waveforms that 
we extracted in this paper with the standard method 
is that we were able to "postprocess" them in order to 
remove an offset at late times. By doing so, we could 
accurately extract the quasinormal frequencies. However, 
we explicitly demonstrated that the large errors in the 
standard method were not due to an overall offset in the 
whole wave. Even after removing the offset "by hand" , 
errors of roughly the same order in the waves remained 
at early and intermediate times in the ringing regime. In 
addition, this postprocessing made use of the fact that 
we knew the qualitative behavior of the exact solution in 
the quasinormal ringing regime. In particular, we knew 
that it had to oscillate around zero, and we also knew 
what the frequencies they were supposed to have. It is 
not clear that one could apply such a postprocessing to 
decrease systematic errors in a more general scenario, 
where the characteristics of the expected waveforms are 
cither completely unknown or not known with so much 
detail. 

Concluding, in this paper we considered weak pertur- 
bations of Schwarzschild black holes, for which — as men- 
tioned above — one can construct the Regge- Wheeler and 
Zcrilli functions (or, equivalcntly, ^4) in an unambiguous 
way everywhere. In a more generic case (say, a collision 
of compact objects) this is not possible, and all gravita- 
tional wave extraction methods are inherently approxi- 
mate at a fixed finite distance. The results of this paper 
suggest that, depending on the accuracy of a given sim- 
ulation, different choices in the extraction procedure at 
a fixed and finite distance may result in relative differ- 
ences in the waveforms that are actually larger than the 
numerical errors of the solution. These differences will in 
general decay with radius, but in a very slow way; typi- 
cally as 1/r (which is, in fact, the decay we found in our 
simulations). For example, in order to decrease the sys- 
tematic errors for an observer at 40 M shown in fig. [7] by, 
say, two orders of magnitude, by just moving the observer 
out and extracting at a single extraction radius, the lat- 
ter would have to be located at ~ 4, 000 M. This means 
that, if similar uncertainties show up in other simulations 
for differing extraction methods, as the results of this pa- 
per suggest (and which can be tested), then decreasing 
those uncertainties by extracting waves at a fixed loca- 
tion and moving the observer further out does not seem 



feasible, and other ideas would have to be explored. 
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Appendix A: VECTOR AND TENSOR SPHERICAL HARMONIC DECOMPOSITION (ODD AND 

EVEN-PARITY SECTORS) 

We discuss now, in some detail, how to compute a multipole decomposition using vector and tensor spherical 
harmonics. A vector field Va defined on the manifold S 2 can be decomposed in multipolcs using even and odd-parity 

basis vectors. Denoting the components in this basis by ft.£ven and h*££ , Va can be written as 

oo I 

Va = E E A&W + hg£>S<£ m \ (Al) 



1=1 m=- 



Here, YY'' m> and S^'^ are the even and odd-parity basis vectors tangent to the sphere, respectively. They are defined 
as 

Y<£' m) = V A Y (e ' m) (A2) 

s (t,m) = gB^ B y(/, TO ) j (A3) 

where Va is the covariant derivative compatible with the unit sphere metric ()ab and cab is the Levi-Civita tensor 

with components Ig^ = sin#. They satisfy the relations VcQab = and Vc^ab = 0. These vectors obey the 
orthogonality relations 

" g AB Yi e ' m) Y^ m,) dn = i(i+l)8 u ,8 mm ,, (A4) 

g AB sf m) S^ m ' ] dn = e(£ + l)6 u >6 mm >, (A5) 

gAB Y ^rn) s (i'm') dri = Q (Ag) 

Here dfl = smd dd d<j) is the area element in polar spherical coordinates and the bar denotes complex conjugation. 
Using the orthogonality property we can find the multipole modes hivcm and . The result is 

a&S? = JaTT) J 9 AB v A n e ' m) dn, (A7) 

*£2° = J WxS£ m) ««■ (A8) 
Using spherical coordinates the components of the even-parity basis are 

Y d {e ' m) = d e Y^ (A9) 

Yj e ' m) = 3^ m \ (A10) 

For the odd parity we get 

s °' m) = ~^te d,pY{£ ' m) (A11) 

Sf m) = waO d e Y^ m \ (A12) 
Expanding the integral with these vector components we obtain that 

= ^ / VeY^ + -A-V,Yj e ^ dn, (A13) 



T 1)J ' sin 2 ., 

For tensors, the idea is the same. If Vab is a tensor field defined on the unit sphere, the multipole decomposition 
takes the form 

oo i 

v ab=Y. E K^ m) g AB Y^ + G& m >Y& m) + hf m) sf™\ (A15) 
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where g A BY {Um) and Y^% m) arc the even-parity tensor basis, whereas S^b^ is the odd-parity tensor basis. We follow 
the Regge-Wheeler notation by using K and G for the even-parity components and for the odd-parity one. The 
tensor basis is defined as 



S 



AB 



VbS 



(A16) 
(A17) 



This definition agrees with Zerilli tensor harmonics up to a factor of 2, as we will see. They obey the orthogonality 
relations 



g Ac g BD Y^ l) Y^'<m = h(£ - 1)(£ + !)(£ + 2)6^ 



g AG g BD S^S A ^'dn = 2£(l -!)(£+!)(£ + 2)8 U >&. 



1 



(A18) 
(A19) 



and integration of the product of two different tensor basis vanishes. With this we can find K, G and li2- The result 
is 



,{l,m) 



^ I V AB g AB Y em dn 



£(£-!)(£ +!)(£ + 2) 



£{£-!)(£ +!)(£ + 2) 
Using spherical coordinates the components of the basis are 

r(t,m 



g Ac g BD V AB Y^dQ 



g AC g BD V AB S^dn 



Y y 
Y ( 

r« 

</). 

S' ( 

si 



2 

2 

1 . 

sin 

2 

1 



2 sin 6 



X 



(£,m) 



2 

-sin<9X^ m \ 
2 



where W^ m ^ and X^ are defined by Zerilli [§3| as 

W (l,m) _ 2 



d 2 e + b(i+l) 



Y { 



X^ m ) = 2d <j> {d e -cot0)Y^ e ' m \ 



(A20) 
(A21) 
(A22) 

(A23) 
(A24) 
(A25) 
(A26) 
(A27) 
(A28) 

(A29) 
(A30) 



Assuming that V AB is a symmetric tensor and abbreviating the normalization constant as L = £(£ — 1)(£ + 1)(£ + 2), 
we expand the integrals to get 



(e,m) 



K 



i / VeeW^ + -2\- (2V H X^ - 
L J sin \ 



,W {t ' m) ) dfl 



2 LJ sin 3 6 



sm( 



(A31) 
(A32) 
(A33) 



The Y lm are normalized with respect to the standard metric g AB on S 2 , an exception being the cases £ = and £ — 1; 
where we choose the normalization such that Y°'° = 1, and J s2 Y 1 ' m Y 1 ' m dfl = 47r/3. 
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